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Abstract 

A holonomic function is a differentiable or generalized function which 
satisfies a holonomic system of linear partial or ordinary differential equa- 
tions with polynomial coefficients. The main purpose of this paper is to 
. present algorithms for computing a holonomic system for the definite in- 

£\J ' tegral of a holonomic function with parameters over a domain defined by 

polynomial inequalities. If the integrand satisfies a holonomic difference- 
O ■ differential system including parameters, then a holonomic difference- 

' differential system for the integral can also be computed. In the algo- 

rithms, holonomic distributions (generalized functions in the sense of L. 
Schwartz) are inevitably involved even if the integrand is a usual function. 
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£J ■ Introduction 

Holonomic systems of linear differential equations, which play a central role in 
the D-module theory, were introduced by Bernstein [2] in the algebraic setting, 
and by Sato et al. (20j in the analytic setting under the name of 'maximally 
overdetermined systems'. We follow the formulation by Bernstein, which would 
be the more adapted to practical applications with computers. Hence, in the 
£C) ' present paper, we mean by a holonomic function a function which satisfies a 

holonomic system of linear differential equations with polynomial coefficients. 
Two equivalent definitions of a holonomic system will be recalled in Section 1. 

Most of the special functions in one variable such as various hypergeomet- 
ric functions and the Bessel function are holonomic by the definition. As an 
important class of holonomic functions in several variables, let us consider 
the multi- valued analytic function u = fi 1 --- with non-zero polynomi- 
als fx, ... , f m in n variables. As a multi- valued analytic function, u is defined 
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on {x £ C" | fx(x) ■ ■ ■ fm(x) 7^ 0} and is holonomic for any complex number 
Xj. We can regard this function also as a distribution (/i)+ • • • (/ m )+ m defined 
on R" in the sense of L. Schwartz if fj are real polynomials and (Ai, . . . , A m ) 
avoids some exceptional set. Such a distribution was introduced by Gel'fand 
and Shilov [7] in some restricted cases. See e.g., [TU] for a theoretical study on 
generalized functions including such a distribution. In particular, substituting 
zeros for the parameters yields 

(fi)°+---(f m )° + = Y(f 1 )---Y(f m ), 
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where Y denotes the Heaviside function, i.e., Y(t) = 1 for t > and Y(t) = 
for t < 0. The Heaviside function will play an essential role in the present paper. 

Algorithmic studies, especially the integration, of holonomic functions were 
pioneered by Almkvist and Zeilberger pQ who introduced what is called the 
creative telescoping method which applies to both difference and differential 
holonomic systems. Then by using Grobner bases in the ring of differential 
operators, Takayama ([23], [24]) introduced two algorithms for integration of 
holonomic functions. These algorithms were generalized to Ore algebras by 
Chyzak [5] (see also 

On the other hand, a precisely Z?-module theoretic algorithm was given in 
[IB] (see also [2T], [H]), which is the 'Fourier transform' of the restriction algo- 
rithm firstly introduced by Oaku 15 and generalized in |19j . Given a holonomic 
system for a function u(x\ 1 . . . , x n ), this algorithm outputs a holonomic system 
for the integral 

v(xi,...,x n - d ) = \ u{xi,...,x n )dx n - d +\---dx n . 

This algorithm might not be efficient enough but has an advantage in the fol- 
lowing two respects: First, the holonomicity of the output is guaranteed in the 
-D-module theory if the input is holonomic. Second, it applies to generalized 
functions as well as to infinitely differentiable functions since it does not involve 
computation with rational functions as coefficients. 

The purpose of the present paper is to apply this D-module theoretic algo- 
rithm to the integral over a domain defined by one or more polynomial inequal- 
ities by using the Heaviside function, generalizing and elaborating a method 
sketched in |17| . More precisely, we are concerned with an integral 

v(xx,...,x n - d ) = / u(x 1 ,...,x n )dx n ^d+i---dx n (1) 

J D(xi,...,x n -d) 

over a domain 

D(xi,...,x n - d ) = {(x n - d+ x,...,x n ) e R d | fi(xi,...,x n ) > 0, . ..,/ m (iri,...,a; n ) > 0} 

with polynomials /i, . . . , f m with real coefficients. 
Let us present two explanatory examples. 

Example 1 Set 

f 1 e x - I 

v(x) = / e xy dy = . 

Jo x 

The integrand u(x,y) = e xy satisfies a holonomic system 

(d x - y)u(x, y) = (d y - x)u(x, y) = 

with d x — d I dx and d y — d/dy. If we apply the integration algorithm to this 
holonomic system, we get an incorrect equation xv(x) = 0, which follows from 
formal integration of the equation (d y — x)e xy = with respect to y; this integral 
does not exist in fact. 

In order to get rid of the boundary conditions at y = 0, 1, we rewrite this 
integral in terms of the Heaviside function Y{t) to get 

/oo 
e x W{y)Y{l-y) dy. 
-OO 
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The new integrand u(x, y) = e xy Y(y)Y(l — y) satisfies a holonomic system 

y(y - l ){dy - x)u{x, y) = (d x - y)u(x, y) = 

in the sense of distribution theory. In fact one has 

y(y - l)(d y - x)(e^Y(y)Y(l - y)) = y(y - l)e^(%) - 5{y - 1)) = 0, 

where 5(y) — d y Y(y) denotes the Dirac delta function. The integration algo- 
rithm applied to this holonomic system outputs an answer 

{xd 2 x -{x- 2)d x - l)v{x) = 0, 

which is correct as is seen by rewriting the operator as follows: 

xd 2 x - (x - 2)d x - l = d x (d x - l)x. 

Of course, one can treat integrals like this one over an interval by the clas- 
sical creative telescoping method taking the boundary conditions into account 
(see e.g, [5]). See also a recent work by [11] for an algorithm to compute in- 
homogeneous differential equations for such an integral. However, it does not 
seem straightforward to apply these methods to more general integrals such as 
the following: 

Example 2 Set 

«(*) = / TT^n: 2 = / (! + - 2 + y^Yit - x 2 - y 2 ) dxdy. 

Note that v(t) is continuous in t £ R and v(t) = holds for t < 0. By using an 
algorithm in |T8], we get generators 

yd x -xd y , (-.t 2 - y 2 + t)d t - s i; 

{{—x 2 + t)d t — s\)d x — yxd t dy — 2xd t , 

(_ x 2 _ y 2 _ l)dy + ( _ 2< _ 2)ydt + (2si + 2s2)y; 

(a; 2 + l)d x + yxdy + (2t + 2)xd t + (-2 Si - 2s 2 )x, 
((-t - l)xd t + Sl x)8 x + {{-t - l)yd t + s iy )d v + (-2t 2 - 2t)8 2 
+((4si + 2s 2 - 2)t + 2si - 2)d t - 2s 2 - 2s 2 si 

of the annihilating ideal of the analytic function (t — x 2 — y 2 ) Sl (1 + x 2 + y 2 ) S2 
with parameters Si,S2- Substituting for s% and —1 for s 2 gives a holonomic 
system for the distribution 

u(x, y, t) = (l + x 2 + y 2 )- l Y{t - x 2 - y 2 ). 

Then applying the integration algorithm to this holonomic system, we obtain a 
differential equation 

{{t 2 +t)d 2 +tdt)v = 0. 

In fact, integration in polar coordinates gives v(t) = irY(t) log(l + t). It is easy 
to verify that v(t) satisfies the above equation as distribution. 
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For the general integral ([lj, we rewrite it in the form 



v(xj, . . ,,x n -d) = / u ( x i, ■ ■ ■ ,x n )Y(f 1 (x 1 ,. . .,x n )) ■ ■ ■Y(f m (xi,. . . ,x„)) dx n ^d+i • • ■ dx n . 

Hence the integration algorithm can be applied if a holonomic system for the 
product uY(fi) . . . Y(f m ) is computed. We shall describe the procedure in detail 
and prove its validity. We remark that the 'indefinite integral' 



is a special case of (Q~]). 

Our algorithm consists of the following two steps: 

1. For a given holonomic system for a function u, compute a holonomic 
system for uY(fi) ■ ■ ■ Y(f m ). 

2. Compute a holonomic system for the integral J Rd uY(fi) ■ ■ ■ Y(f m ) dx n -d+i ■ ■ ■ dx, 

For the first step, we begin with an algorithm to compute a holonomic system 
for Y(fi) ■ ■ ■ Y(f m ). Then the tensor product computation for Z?-modules gives 
an answer to the first step. For the case where u is a complex power, or an 
exponential of a polynomial, the tensor product computation is unnecessary 
as in the previous examples. Even for general u, we shall give an alternative 
method which avoids the tensor product computation. The second step can be 
done by the integration algorithm for D-modules. 

In many practical examples, the integrand can have auxiliary parameters 
other than x\,...,x n above but cannot be regarded as a holonomic function 
including the parameters. For example, consider the integral 



with parameters s±, . . . , s m . The integration algorithm cannot be applied di- 
rectly unless we specify the values of the parameters si, . . . , s rn explicitly. How- 
ever, it is often the with the example above, that the integrand satisfies 
a holonomic system of difference-differential equations including the parameters. 
We also give an algorithm for computing a holonomic difference-differential sys- 
tem for such an integral. 

We have implemented the algorithms in the present paper by using a com- 
puter algebra system Risa/Asir [T^]. In particular, we make use of a library 
file nkjrestriction by H. Nakayama (see for computing restriction and 
integration. 

1 Differential operators and holonomic systems 

Let us denote by D n the ring of differential operators on the variables x — 
(xi, . . . ,x n ) with polynomial coefficients. An element P of D n is written in a 
finite sum 
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where a — (ai, . . . , a n ), (3 = . . . , (3 n ) G N™ are vectors of nonnegative 
integers with N = {0,1,2,...}, x a = x" 1 = d{ 1 • • • dfr with the 

derivations di ~ d/dxi (i = 1, . . . , n), and a Qjj a are complex numbers. 

Given Pi, . . . ,P r G D n , we associate the left ideal / — D n P\ + • • • + D n P r 
generated by P\ , . . . , P r with a system of linear differential equations 

Piu = ■■■ = P r u = (3) 

for an unknown function u. This enables us to work with a left ideal of D n 
instead of each system of linear differential equations. Here we suppose that the 
unknown function u belongs to a 'function space' T which is a left D„-module. 
Examples of such T are the set C°°(U) of C°° functions on an open subset U of 
R", the set 0(U) of possibly multi- valued analytic functions on an open subset 
U of C n , the set V'(U) of the Schwartz distributions on an open subset U of 
M", and the set <S"(R ra ) of tempered distributions. 
A weight vector for D n is an integer vector 

w= (wi,...,w n ;w n+1 ,---W2n) € 

with the conditions Wi + w n +i > for i = 1, . . . , n, which are necessary in view 
of the commutation relation diXi — Xidi + 1 in D n . For a nonzero differential 
operator P of the form ([2]), we define its w- order to be 

ord w (P) := max{(w, (a,/3)) = wiaiH hm„a„+w„ + i/3iH Vw 2n Pn \ a a ,0 0}- 

We set ordu,(0) := — oo. A weight vector w induces the w -filtration 

F?(D n ):={PeD n \ord w (P)<k} (k G Z) 

on the ring D n . Following Bernstein [2], let us define the notion of holonomic 
system by using the weight vector (1,1) = (1, . . . , ; 1, . . . , 1) <E l? n . Let M be 
a left -D ra -module and {F/ c (M)}/ ce z be a good (1, l)-filtration. This means the 
following properties: 

1. every Fk(M) is a finite dimensional vector space over C; 

2. F k {M) C F k+ i(M) for all k e Z; 

3. |jF fe (M) = M; 

fcez 

4. F^iD^FkiM) C ^ +fc (M) for all i, k e Z; 

5. there exists fci g Z such that Fk(M) = for < fci; 

6. there exists fc 2 G Z such that \D n )F k {M) = F i+k (M) for fc > fc 2 . 

Then there exists a polynomial in fc such that dime F k {M) = p(k) for sufficiently 
large k. The degree of p(k) does not depend on the choice of a good (1, 1)- 
filtration of M and is called the dimension of the module M, which we denote 
by d(M). It was proved by Bernstein [5] that d(M) > n if M ^ 0. Following 
Bernstein let us adopt the following 
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Definition 1 A finitely generated left D„-module M is called a holonomic sys- 
tem if d(M) < n. We also call a left ideal I of D n to be a holonomic ideal, by 
abuse of terminology, if the left D„-module D n /I is holonomic. 

Note that d(M) < n is equivalent to d(M) = n or M = in view of the 
Bernstein inequality stated above. The dimension d(M) can be computed as 
the degree of the Hilbert function from a Grobner base with respect to a term 
order which is compatible with the total degree. 

Let us recall another characterization of a holonomic system essentially 
given by Sato et al. [20]. For this we use the weight vector w = (0,1) = 
(0, . . . , 0, 1, . . . , 1). Let P be a nonzero differential operator written in the form 
@ and set m :— ord/ ,i)(-P)- Then the principal symbol of P is the polynomial 
defined by 

a(p)( X ,o= E E a ^ a ^ 

|/9|=m a 

where £ = (£i,...,£ n ) are the commutative variables corresponding to the 
derivations d = (di, . . . , d n ). 

Set M :— D n /I with a left ideal / of D n . The characteristic variety of M is 
defined to be the algebraic set 

Char(M) := {(x, £) € C 2n | a(P)(a!, 0=0 for any P € I\ {0}} 

of C 2 ™. It was proved in [20] that the dimension of (every irreducible component 
of) Char(M) is not less than n unless M — 0. Especially, M is holonomic if 
and only if the dimension of the characteristic variety is n or else M = 0. 
The equivalence of these two definitions of holonomic system is proved, e.g. in 
Chapter 3 of [3]. The characteristic variety can be computed via a Grobner 
base with respect to a term order which is compatible with the (0, l)-weight 

(cf. na). 

2 Holonomic distributions 

First let us recall the definition of distributions due to [22 . 

Definition 2 Let Cq°(U) be the set of the C°° functions on an open set U of 
K™ with compact support. A distribution u on U is a linear functional 

u:C^{U) 3^(«,^)eC 

such that lim j oo (u, (fij) = holds for a sequence {ifij} of C^°(U) if there is a 
compact set K contained in U such that <pj are zero on U \ K and 

lim sup \d a (fij(x)\ = for any aeN". 

i^oo X<£U 

The set of the distributions on U is denoted by V(U). The derivative d^u of u 
with respect to Xk is defined by 

(d k u, ip) = -(u, d k (fi) for any ip e Cfi°(U). 

For a C°° function a on U, the product au is defined by 

(au, ip) — (u, aip) for any p £ Cq°(U). 
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In particular, by these actions of the derivations and the polynomial multipli- 
cations, D'(U) has a natural structure of left D„-module. 

Definition 3 Let u be a C°° function or a distribution defined on an open 
subset U of K™. Then we call u a holonomic function or a holonomic distribution 
on U if u satisfies a holonomic system. In other words, u is holonomic if and 
only if its annihilator 

Axm Dn u := {P G D n \ Pu = on U} 

is a holonomic ideal. 

Let fx, . . . , f m be polynomials with real coefficients. We assume that the 
set {x G 1" | fi(x) > (i = l,...,m)} is not empty. Then the distribution 
u = (A)^ 1 • • • (/ m )+ m on R n is defined to be 

(u,y)=( A(x) Al ■■■ f m {x) x ™ V (x)dx 

•>/i>0,...,/ m >0 

for (p G C^°(M™) if Re Ai > for each i. Moreover, u, that is, (u,(f) for any 
tp G Co°(R n ), is holomorphic in (Ai, . . . , A m ) on the domain 

fl+ := {(A 1; . . . , A m ) g C m | Re A 4 > (i = I, . . . , m)} 

and is continuous in (Ai, . . . , A m ) on the closure of f2+. 

In order to deduce a holonomic system for the distribution (/i)+ ■ • • (/ m ) + m , 
let us first work with the 'formal' function /-j* 1 • • • /^ m . More precisely, setting 
F = fl ' ' ' /m j we consider a free module 

£:=C[a;,,s,^- 1 ]/r •••/£" 

over C[x, s, -F^ 1 ], which has also a natural structure of left _D„[s]-module induced 
from the formal derivation 

wv ■ ■ • m = ^/r • • • /to- + «E^H//Vr • • • /*» 

1 j=i 1 

with a G C[ar, s,.F ]. Then the annihilator ideal 

Ann Cn[s] /f •••/'•» = {P(.s) G D n [«] | P(s)/f •••/*»= in £} 

can be computed precisely with an algorithm of [T5], or of [I] if the coefficients 
of /i, . . . , / m are contained in a computable subfield of C. 
Define a Z?„[s]-submodule A/" of C by 

:= D n [a]f? ■ ■ ■ f£ =s DnM/Aim^w/f 1 ■ ■ • /£" . 

Let us specialize the parameters s. For A = (Ai, . . . , A m ) G C' m , set 

7V(A) := AA/(( Sl - Ai)AT + • ■ • + (.s m - X m )J\f). 

Then A/"(A) has a natural structure of left D„-module which is induced from the 
left D n [s]-module structure of N by the identification 

D n ~ D„[s]/((si - Ai)D n [s] + • • • + (s m - A m )D n [s]). 
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Proposition 1 Af(\) is a holonomic D n -module for any AeC". 

Proof: We use a filtration 

k 

Fi(D n [s]) := {P(s) = J2 P J sj I P J e fl n, i + ord (lil) (P,) <i,kN} (i e Z) 

on Z? n [s]. Let us define a filtration 

F k (C):={±f? \ aeC[x,s],dega<dk} (k e N) 

on £, where deg denotes the total degree in (x, s) and = degF + f . Then it 
is easy to see that 

F k (Af) := F k (D n [s])f? ■ ■ ■ /*» C F fe (£) 
holds for fc > 0, which implies an inequality 



dim c F fc (AA) < dim c F fc (£) 



dk + n + m 
n + m 



There exists a polynomial p(k) in k such that dime F k (J\f) = p(k) holds for 
integers k large enough since {F k (J\f)} is a good filtration on N. From the 
above inequality, it follows that the degree of p(k) is at most n + m. 
Set Af'(X m ) '■= N / (s rn — X m )Af and define a good filtration on it by 

F k (N'{\ m )) = F k (D n [ s })f s /(F k (D n [s})f s n (s m - X m )Af) (k e N) 

with f s := • • • Since (s m -A m )F fc _i(L>„[s])/ s is contained in F fe (D„[s])/ s n 
(s m — X m )Af and the homomorphism s TO — A m : Af — > .A/" is injective in view of 
the definition of A/", we have 

dim c F fe (Jv"(A m )) = dim c F k (D n [s})f s - dim c {F k {D n [s])f s n (s m - A m )W) 
< dim c F fe (£)„[ S ])/ s - dime Jlfc-i^nW)/ 8 
= p(A:)-p(fc-l) 

for k sufficiently large. It follows that dime F k (J\f'(X m )) coincides with a poly- 
nomial in k of degree < n + m — 1 for sufficiently large k. Proceeding in the 
same way we can show that dime F k (J\f(X)) is a polynomial in k of degree < n 
for k sufficiently large with 

F k (Af(X)) := F k (Af)/(F k (Af) n ((si - \i)Af + ■ ■ ■ + (s m - X m )Af)). 

Note that we have 

F k (Af(X))=F^ 1 \D n )[f^...f^} 

with [/f 1 • • • /4 m ] being the modulo class of f^ 1 ■ ■ ■ in 7V(A). This implies 
that AT(X) is a holonomic £>„-modulc. □ 

Now let us make explicit the relation between the algebraic module Af(X) 
and the distribution (/i)^ 1 • • • (/ m )+ m . 
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Lemma 1 Let P(s) = P(si, . . . , s m ) be an element of D n [s\. Then 

P(Xl,...,X m )(fl)+ 1 ---(fm) > r =0 
holds m V'(R n ) for any (Ai, ...,X m ) € fl + if and only if P(s)ff • • • /*» = 

Proof: It is easy to see that 

m f)f 

holds as distribution and the right-hand side is a locally integrable function if 
Re Xj > 1 for j = 1, . . . ,m, in accordance with the action of di on C. Let 
P{s) £ D n [s]. There are a(x, s) £ <C[x, s] and a positive integer k such that 



holds in £. Moreover, the right-hand side of this equality vanishes in C if and 
only if a(x, s) — 0. From the observation above, we know that 

P(X U . . . , X m )(h)^ ■ ■ ■ {f m )\- = a(x, A)(/!)^- fc ■ • • {f m ) x +- k 

holds if Re Xi > k {i — l,...,m). Since the distribution on the right-hand 
side is a locally integrable function, it vanishes if and only if a(x, X) — 0. The 
conclusion follows from the uniqueness of analytic continuation. □ 

Proposition 2 If the distribution (/i)! 1 • • • (/m)+ m is well-defined and holo- 
morphic in s = (s\, . . . , s m ) on a neighborhood of X = (Ai, . . . , A TO ) 6 C m , then 
there is a surjective D n -homomorphism 

N(X) — > ■ ■ • {fm)+ m 



which sends the residue class [/j 51 • • • <E A/"(A) to (/i)^ 1 ■ • • (/ m )+ m - 
Proof: Let us define a surjective 23 n -homomorphism $ : Af(X) — > D„(/i)^ 

by 

HP[f! 1 ---f% n ]) = P(A)\ 1 ---(fm) > ? 

for P e D n . If P[/f • • • = in Af(A), then there exist Pi(s), . . . , P m (s) € 
D n [s] and Q[s] e Annfl^^/j 1 • • • such that 

P = ( Sl - Ai)Px(s) + • • • + (s m - X m )P m (s) + Q(s). 

Then in view of Lemma [1] and the uniqueness of analytic continuation, we have 

P(fl)^ ■ ■ ■ = QW(fl)+ ■ ■ ■ (fm)+ m = 0. 

Hence $ is well-dehned as a homomorphism of left D„-modules. □ 

It is known that there exist a non-zero polynomial b(s) in s — (s%, . . . , s m ) 
and an operator P(s) € -D n [s] such that 

V 

P(s)/r +1 • • • f s n T +1 = 6(s)A Sl • • • /*», = n( c « s i + ' ' -+Cims m + Ci) (4) 
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with positive integers and positive rational numbers Cj. This was proved 
by Kashiwara [3] for m = 1 and by Gyoja [5] for general m. An algorithm to 
compute the ideal consisting of such b(s), which is called the Bernstein- Sato 
ideal, was given in [T5]. For the case m = 1, an algorithm to compute such a 
functional equation was given in [Tl] . See also [TB] for an algorithm to compute 
a minimal one in the sense that the differential operator on the left-hand side 
has minimal degree with respect to d x and s. 
In view of Lemma [TJ we obtain from (Tj[]) 

ft(Ai, . . . , A m )((/ 1 )^ 1 • • • (f m )l m ,<f) = (P(Ai, . . . , A m )(/i)^ 1+1 • • • (/ m )^ +1 ,^) 

= ((/i)i 1+1 ---(/ m )t m+1 , t J P(A 1 ,...,A m V) 

for </? e C£°(R n ) with the formal adjoint operator ( P(s) € D„[s] of P(s). By 
using the functional equation (g]) repeatedly, we know that • ■ ■ (/ m )+ m can 

be continued to a distribution in x which is meromorphic in (Ai, . . . , A m ) on the 
whole C m . More precisely, (/i)^ 1 • • • (/ m )+ m is holomorphic in (Ai, . . . , A m ) on 

n(/i, ■ • • , f m ) ■= {(Ai, . . . , A m ) £ C m I 6(Ai + fc, . . . , \ m +k) + for any k £ N}. 

Lemma 2 The set 0(/i, . . . , / m ) contains the closure of Q+. 

Proof: Assume that A' - 1 = (A^ , . . . , A^') belongs to the closure of i.c, 
Re A[ ' > for any i. Then we have 

Re (ca(A^ 0) +&) + ••• + Cim(A^ + fc) + c,) > c, > 

for any k £ N. This implies 6(A^ 0) + fc, . . . , \$ + k) ^ 0. Hence A^ belongs to 
fi(/i,...,/ m ). □ 

Propositions [1] and [21 provide us with an algorithm to compute a holonomic 
system for (/i)^ 1 • • • (/ m )+™ , and hence for Y(fi) ■ ■ ■ Y(f m ) by virtue of Lemma 
® 

Algorithm 1 (a holonomic system for (fi)+ ■ • • (/m)^™) 

Input: fx, . . . , f m £ R[x] and A = (Ai, . . . , A m ) £ Q(/i, . . . , f m ). 
Output: A set G of generators of a holonomic ideal contained in 
AnnoMi)^ ■ ■ ■ (fm) X + m ■ 

1. Compute a set G\ of generators of the annihilator Anii^fji/J 1 ■ • • / r ^ m by 
an algorithm of [18] . or of [1]. Following the former, let J be the left ideal of 
the ring of differential operators on the variables x±, . . . , x n , t\ , . . . , t m which is 
generated by 

m df 

Then compute the intersection J D D n [s] with the identification sj — —dt tj 
(J = 1, . . . , m) by using Algorithm [7] in Section 5. 

2. Set G := {P(A) | P(s) e Gi}. 
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We do not know if the output of this algorithm coincides with the annihilator 
ideal of (/i)+ • • • ( f m )+ m ■ As a trivial counter example, take / = x 2 with a 
variable x. Then d x Y(f) = d x l — holds. On the other hand, the annihilating 
ideal of f s = x 2s is generated by xd x — 2s. Hence d x cannot be obtained by 
substitution s = from an annihilator of f s . 

Example 3 Set / = x 3 — y 2 . From the functional equation 

(s* + \< - ^) ^ = < A + •> ( A + 1) ( A + 1) 

we know that the distribution u := is holomorphic in A belonging to the set 

n(/) = {A€C| A^-l-fc, -^-fc (A = 0,1,2,...)}. 

A holonomic system for u is given by 

(2yd x + 3x 2 d y )u = (2xd x + 3yd y - 6X)u = 0. 

In particular, the Heaviside function Y(f) satishes a holonomic system 

(2yd x + 3x 2 d v )Y(f) = (2xd x + 3yd y )Y(f) = 0. 

Differentiation of the distribution (/i)+ • • • (/ m )+ m with respect to the com- 
plex parameters Ai, . . . , A m yields a distribution 

(/^■■■(/^(log/ir--- (log 

which is holomorphic in A on . . . , f m ) f° r arbitrary non-negative integers 

i/i, . . . , z^ m . A holonomic system for this distribution can be computed through 
differentiation of the annihilator of f* 1 • • • with respect to the parameters 

^1 5 • • * ; • 

For the sake of simplicity, let us describe an algorithm in case m = 1. The 
general case m > 2 is quite similar. 

Algorithm 2 (a holonomic system for /+(log/) fe ) 

Input: / e M[x], fc e N and A € 

Output: a set G of generators of a holonomic ideal contained in 
Ann D „/ A (log/) fe . 

1. Compute a set Gi of generators of the annihilator Ann£> n [ s ]/ S . 

2. Let ei = (1, 0, . . . , 0), • • • , et+i = (0, . . . , 0, 1) be the canonical base of Z fc+1 . 
For each P(s) € Gi and an integer j with < j < k, set 

i=0 v 7 

3. Let A be the submodulc of Z?„[s] fc+1 generated by the set 

{P^\s) | P(s) <E Gi, < j < k}. Compute a set G2 of generators of the ideal 

I := {Q(«) G D n [s] I Q( S )e fc+1 e A} 

by eliminating ei, . . . , e& via a Grobner base of A with respect to a 'position 
over term' ordering. 

4. Set G := {(2(A) | Q(s) e G 2 }. 
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Proposition 3 The ideal I of the preceding algorithm is holonomic and anni- 
hilates the distribution /+(log/) fc for any A G ^(/)- 

Proof: Let P(s) belong to G\. Differentiating the equation P(s)f s = with 
respect to s, one gets 

i=0 v 7 

Hence P^(A) annihilates the vector log /,..., /+(log f) k ) of distribu- 

tions for A G It follows that each element of G annihilates /+(log/) fe . 

In order to prove the holonomicity, let N(X) be the submodule of (D n ) k+1 
generated by {P(A) | P(s) G N}. Then it suffices to prove that (D n ) k+1 /N(X) 
is holonomic since D n /I is its submodule. We regard (D n ) J+1 as the submodule 
of (D n ) k+1 generated by ei, . . . , e J+ i if < j < k. Set 

Af(XJ) := (D n y +1 /(N(X) n (£)„y +1 ). 

Then we have an increasing sequence 

W(A, 0) c Af(A, 1) c • • • c jV(A, k) = {D n ) k+1 /N(X) 

of left Z?„-modules. We have only to show that 7V(A, 0) and Af(X, j)/N(X, j — 1) 
are holonomic for 1 < j < fc. 

Note that Af(\), which annihilates fi_, is isomorphic to D n /I(X) with 

/(A) := {P(A) | P(s) G Ann Dn [ s ]f s }. 

From the dehnition we have I(X) C iV(A) n D n . Together with Proposition 
[T] this implies that Af(\, 0) is holonomic. Let ij be the homomorphism of D n 
to (D n y +1 which sends P G D n to Pe J+ i. It is easy to see that I(X)ej+i is 
contained in (N(X) fl (D„) J+1 ) + (D n y . Hence t 3 induces a surjective homomor- 
phism of the holonomic module Af(X) to Af(X,j)/AT(X, j — 1), which is hence 
also holonomic. This completes the proof. □ 

Example 4 Set f = x 3 — y 2 . Then Algorithm [5] outputs three annihilators 
2yd x + 3x 2 d y , 

Ax 2 d 2 x + (Uyxdy + (-24A + 4)x)d x + 9y 2 d 2 + (-36A + 9)yd y + 36A 2 , 

- 8yd% + (36t/x9| - 72Xxd y )d x + 27y 2 d 3 + (-108A + 8l)yd 2 + (108A 2 - 108A + 3)5, 

of /^log/. 

Finally let us consider the delta function 5(f) for a non-singular real polyno- 
mial / G R[x] in n variables. As distribution, 6(f) is defined to be the integral 




tpu) 



for tp G Cg°(R n ), where w is the volume element of the hypersurface / = such 
that df A uj — dx\ A • • • A dx n . The following should be well-known: 
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Proposition 4 Let f e M.[x] be non-singular in C n ; i.e., the variety defined by 

{i6 ch /«-£«-■■ •-£«-<>) 

&e empty. Then 5(f) satisfies a holonomic system 

Proof: The characteristic variety is contained in the n-dimensional algebraic set 



G C 2 " I /(*) =0,6 = c^-(i = l »),cg 



Hence the system above is holonomic. □ 

3 Operations on holonomic functions 

Let us assume that u and v are holonomic functions or distributions. Then the 
following functions or distributions are holonomic under the condition that they 
are 'well-defined'. Moreover, if a holonomic ideal for u (and also one for v if 
relevant) is explicitly given, then there exist algorithms to compute a holonomic 
ideal which each of the following functions satisfies. 

1. Pu with P e D n . 

2. The sum u + v. 

3. The restriction u\y of u to an affine subspace Y of K™ if it is well-defined 
as in the case where u is smooth. 

4. The product uv if it is well-defined as in the case where u is C°° and v is 
a distribution. 

5. The definite integral / u(x\, . . . , x„_ d , x n -d+i, ■ ■ ■ , x n ) dx n - d+1 ■ ■ ■ dx n 

JR d 

with parameters if it is well-defined as in the case where u has a compact 
support with respect to the integration variables. 

Let us explain briefly how to compute holonomic ideals for functions above. 
The algorithms which will be presented in this section are more or less well- 
known at present. Nevertheless, we would like to recall them in order to verify 
that they certainly apply to holonomic distributions as well as for the reader's 
convenience. 

Let / and J be holonomic ideals for u and v respectively. First, a holonomic 
ideal for Pu can be computed as an ideal quotient / : P by using Grobner 
bases in the same way as in the polynomial ring. The left _D„-homomorphism 
of D n to itself which sends Q € D n to QP induces an injective homomorphism 
D n /(I : P) -> D n /I. Hence D n /(I : P) is holonomic. 

Second, a holonomic ideal for u + v can be computed as an ideal intersection 
I (1 J. The left _D„-homomorphism of D n to (D n ) 2 which sends Q <G D n to 
(Q, — Q) induces an injective homomorphism 

D n /(Ir\J)^(D/I)®(D/J). 
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This implies that D n / (I n J) is holonomic. 

The restriction algorithm was given in |15j for one codimensional case and 
in Q]5] for the general case. The restriction algorithm is translated to the in- 
tegration algorithm through the algebraic Fourier transformation of the ring of 
differential operators. Hence let us describe the integration algorithm instead. 

Let u(x,t) be a holonomic distribution defined on (7 x R d in the variables 
(x, t) — (x%, . . . , x n , ti, . . . , td), where U is an open set of W 1 . Then the definite 
integral of u{x, t) with respect to t = (t\, ... ,td) is defined as the distribution 

C^(U) 3 tp(x) .— ► ( f u(x,t)dt!-"dt d , if) = (u(x,t),tp(x)l{t)), 

where l(t) denotes the identity function which takes only the value 1 for all 
t. This definite integral is well-defined if u(x, t) is a distribution with proper 
support with respect to t, i.e., for each xq £ U, there exists a neighborhood V 
of xq and a compact set K of R d such that u(x, t) vanishes on V x (R d \ K), or 
else if u(x,t) is C°° in (x,t) and rapidly decreasing with respect to t, i.e., for 
any Xq and a, (3 £ N d , there exists a neighborhood V of £o such that t a d^u{x, t) 
is bounded on V x W 1 . 

Let J be a left ideal of D n+ d which is contained in Ami£, n+d u(x, t). The 
D-module theoretic integration of the module D n+ d/I is defined to be the left 
Djj-module 

J D n+d /Idt := D n+ d/{d tl D n+d + ■■■+ d td D n+d + I). 

Let w = (wi, . . . , W2n+2d) € Z 2 (™ +d ^ be the weight vector for D n+ d such that 
w n+ i = 1 (the weight for ti) and u>2n+d+i = — 1 (the weight for <9 ti ) for i — 
l,...,d and that other components of w are zero. If I is a holonomic ideal, 
then there exists a nonzero univariate polynomial b(s) (the 6-function of / with 
respect to w) of the minimum degree such that 

b(-d tl h d td t d ) + PeI 

with some P £ F^ 1 (D n+ d). Let k\ be the maximum integral root of b(s) = 0. 
(Set ki = —1 if there is none.) Then the D ra -module J D n+ d/I dt is generated 
by the residue classes of t a with a £ N d such that \a\ < k%. Set 

N := {(P a ) a £ D n \ P a t a £d tl D n+ d + ---+d td D n+d + I}. 

The 6-function b(s) and a set of generators of iV can be computed by a Grobner 
base of I with respect to the weight vector w. Then we have an isomorphism 
fD n+d /Idt~ (© |a 

|<fei Dn \ /N, which is a holonomic _D„-module. 

Definition 4 The integration ideal of a left ideal / of D n+ d is defined to be the 
left ideal 

N := {P £ D n \ P £ d tl D n+d + ■■■ + dt d D n+d + 1} 

of D n . 

Once a set of generators of N is obtained, the integration ideal Nq can be 
computed by elimination in the free module ©i^io^-Dri- 
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Theorem 1 Suppose that the definite integral v(x) :— J Rd u(x,t)dt\ ■ ■ ■ dtd is 
well-defined as a distribution on an open set U o/K™ and I is a holonomic ideal 
of D n+C i annihilating u(x,t). Then the integration ideal No of I is holonomic 
and annihilates v(x). 

Proof: Let P = P(x, d x ) £ D n belong to Nq. Then there exist Qi — Qi(x, t, d x ,dt) £ 
D n+ d and R E I such that 

P = d tl Q 1 + --- + d td Q d + R. 

Thus for any <p{x) £ Cg°(U) one has 

(P / u(x,t)dh---dt d , <p) = (Pu(x,t), <p(x)l(t)) 

d 

= ^2(d ti Qiu(x,t), <p(x)l(t)) 
i=i 

d 

= -J2{u(x,t), t Q i d t M*)Ht)))=o, 

i=l 

where t P denotes the formal adjoint of P. 

The Z)„-module D u /Nq is holonomic since it can be regarded as a submodule 
of the holonomic module ©| Q |< fcl D n /N, the proof of the holonomicity of which 
is given, e.g, in [3]. □ 

The integration algorithm is summarized as follows. See [TB] or for the 
correctness proof. 

Algorithm 3 (integration ideal) 

Input: A set Go of generators of a holonomic ideal I of D n+ d annihilating 
u(x, t). 

Output: A set G of generators of a holonomic ideal No of D n annihilating 

v{x) := J Rd u(x, t) dtf- dt d . 

1. Compute a Grobner base G\ of / with respect to a monomial order which is 
compatible with the weight vector w — (0, . . . , 0, 1, . . . , 1; 0, . . . , 0, —1, .... —1) 
for the variables (x, t,d x ,d t ). 

2. Compute a 6-function of / with respect to w, which is a univariate 
polynomial b(s) of the minimum degree such that 6(— dt t tx — ■ ■ ■ — dt d t d ) + P 
belongs to I with some P € F^ 1 (D n+ d). This amounts to computing a 
generator of the intersection ideal 

C[9 tl t a + • • • + d td t d ] n *r(J)/*T-i(J). 

fcez 

3. Let k\ be the maximum integral root of b(s) = if any; if there is none or 
else fci < 0, then set G := {1} and quit. 

4. For P 6 Gi and a 6 N d such that ord„,(P) + \a\ < ki, one has an 
expression of the form 

d 

j=l |/3|<fei 



15 



with Qj G D n+d and R fj G D n . Set x(i"P) := (R/3)\p\<ki- Lct N be the 
submodule of ©i^^^ 7>„ generated by 

{ X (t a P) \P€Gi, OYd w (P) + \a\ < h}. 

5. Compute a set G of generators of N n 7>„ via a Grobner base of N with an 
appropriate 'position over term' ordering, where D n is identified with the 
submodule (B\ p \< D n of ©i^^ D n . 

Finally let us describe an algorithm for the product. Let u(x) and v(x) be 
holonomic distributions defined on K™ and assume that the product (Pu(x))(Qv(x)) 
is well-defined as a distribution for any operators P,Q e D n and satisfies the 
following 'associativity' condition: 

{a{x)Pu){Qv) = (Pu){a{x)Qv) (VP,Q € D n , Va(x) € C[a;]). (5) 

This is certainly the case, for example, when either u(x) or v(x) is C°°. 
Let J and J be left ideals of D n such that 

7 C Ami£> n u, J C Aniifl n i! 

and set M := D n /I and TV := D n /J. Let w and w be the modulo classes of 
1 G D n in M and TV respectively. Then there is a natural bilinear map 

^ : M x N 3 (Pu, Qv) i — ► (Pu)(Qu) G £> n (m>) C P'( R ") 

with the property ^(aPu,Qv) = ^(Pu,aQv) for any a G C[x]. Hence by the 
universal property of tensor product, there is a D„-homomorphism 

$ : M ®c[ x ] N -> D„(w) 

such that $((Pm) ® (Qu)) = V(Pu, Qv) = (Pu)(Qv) for any P,Q £ D n . More- 
over, M (g>c[x] N has a natural structure of left 7>„-module and $ gives a ho- 
momorphism of 7)„-module. This means that P(u <£>v) — in M ® TV implies 
P(m>) = in P'(E n ). It is also well-known that M ® C [x] N is holonomic if M 
and TV are holonomic. (See e.g., [3])- In conclusion we have 

Proposition 5 Lei it, v be distributions which satisfy Lei 7 and J 6e ZioZo- 
nomic left ideals of D n which annihilate u and v respectively. Then the left 
ideal 

Ann_D n (u ® 5) := {P G D n \ P(u (g> v) = in M (g> C [ x] N} 
is holonomic and annihilates the distribution uv. 

The 'algebraic' annihilator Ann^ (u ® v) can be computed as follows: First 
let 1(g) J be the left ideal of 7?2„ which are generated by {P(x,d x ) | P G 7} 
and {Q(y,d y ) | Q G J}. Here we use the notation d x = (d/dxi, d/dx n ) 
and d y — (d/dyi, . . . ,d/dy n ). We change the variables (x,y) to (x,z) by the 
substitution x — x, z — y — x. Then the derivations are transformed by 

d Xi = d Xi -d Zi , d Vi =d Zz (i = 1, . . . , n). 
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We regard I<E)J as a left ideal of the ring of differential operators D X:Z with 
respect to the variables (x, z). Then the restriction of D X>Z /I®J to z = is 
defined to be 

D X<Z /{I®J + Zl D x . z + ■■■+ z n D x . z ) 

as left £>„-module, where D n stands for the ring of differential operators with 
respect to the variables x = {x\, . . . , x n ). One has an isomorphism 

M ® C [x] N ~ D XtZ /(I®J + z x D x . z + ■■■ + z n D x . z ). 

Set 

n 

i<e>J\a ■= {P(x, d x ) e D n I P(x, d x )~Y, ZiQi e I® J (3Qi, . . . , Q„ e 

i=l 

This ideal can be computed by the restriction algorithm and coincides with 
Ann£) n (u ® v) in M ®c[x] N. See [H] for details. 

For practical computations of restriction, integration and tensor product, we 
have made use of a library file nk_restriction.rr written by H. Nakayama for 
a computer algebra system Risa/Asir [T^] . 

Example 5 Let us consider the integral 

v(t) := [ fe 9 S(t ~x 2 -y 2 - z 2 ) dxdydz = Y(t) [ fe 9 uj t 
Jr 3 Js? 

for polynomials /, g, where u t denotes the canonical volume element of the 
sphere S 2 := {(x,y,z) e I 3 | i 2 |j/ 2 + z 2 = t}. The annihilator ideal I of 
8(t — x 2 — y 2 — z 2 ) is generated by t — x 2 — y 2 — z 2 , d x + 2xd t , d y + 2yd t , 
d z + 2zdt- Then the annihilator of fS(t — x 2 — y 2 — z 2 ) is given by the ideal 
quotient I : f. Finally, the annihilator of fe 9 8{t — x 2 — y 2 — z 2 ) is given by 
the ideal e 9 (I : f)e~ 9 . For example, if / = 1 and g — x — y 2 — z 2 , we get a 
differential equation 

(Udf + (At + 6)df + hd t - l)u(t) = 0, 

which has a regular singularity at t — with exponents 0, i, 1. 

On the other hand, if / = 2x 2 — y 2 — z 2 and g = 0, we get v(t) = 0. 

Example 6 The integral 

v (t) := f S(t -x 3 + y 2 )e- x2 ~ y2 dxdy 
Jw 2 

satisfies a differential equation 

(108t 2 <9 t 5 + (-216i 2 + 648t)d? + (108i 2 - 972t + 627)9 t 3 

+ (356i - 606)9 t 2 + (-64* + 108)<9 t + 32t - 48>(i) = 0, 

which has a regular singularity at t = with exponents 0, 1, 2, — ^, |. 
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4 Definite integrals with the Heaviside function 

Let u be a holonomic function on W 1 . Let /i, . . . , f m be nonzero polynomials 
in x = {x\, . . . , x n ) with real coefficients. Setting 

D(xi,...,x n - d ) = {(x n - d+ i,...,x n ) € R d | fi(x) > 0,...,f m (x) > 0}, 

let us consider the definite integral 

v(xx,...,x n -d) = / u(x) dx n - d+1 ■ ■ ■ dx n 

J D(xi,...,x n -d) 

= / Y(fi)---Y(f rn )u(x)dx n -d+i---dx„. 

JR™ 

We suppose that the set {fj(x) > | j = l,...,m} is non-empty and the 
integral is well-defined if (xi, . . . , x n -d) belongs to an open set U of R"~ d . This 
is the case, for example, when u is a C°° function on an open set W of K n such 
that 

{(xi,...,x n - d )} x D(xi,...,x n -d,) C W 

and D(xi, . . . , x„_d) is compact if {x\, . . . , x n ~d) G f ■ Then a holonomic ideal 
for this integral can be computed by combining the algorithms explained in the 
preceding two sections. In conclusion, we have the following algorithm. 

Algorithm 4 (a holonomic ideal for a definite integral) 

Input: a holonomic ideal J annihilating the function u. 
Output: a holonomic ideal J annihilating the definite integral v. 

1. Compute a holonomic ideal annihilating Y{f\) ■ ■ ■ Y(f m ) by AlgorithmQ"] 

2. Compute a holonomic ideal / annihilating the product uY(fi) ■ ■ ■ Y(f m ) by 
using the tensor product computation. 

3. Compute the integration ideal of / with respect to x n -d+i, ■ ■ ■ , x n by 
Algorithm [3] 

In the step 2 of this algorithm, if u is of the form u = (.91)+ ■ ■ • {dp)'+e h with 
polynomials gi, . . . , g p , h, then we can compute first a set G which generates a 
holonomic annihilating ideal I of 

(/i)° + ---(/ m )° + (3i)i 1 ---(ffp)+ P 

by using Algorithm [TJ Then an annihilating ideal for uY{f\) ■ ■ -Y(f p ) can be 
computed as an ideal e h Ie~ h of D n . A set of generators of this ideal is obtained 
by substituting di — dh/di for di (i = 1, . . . , n) in each element of G. 

Example 7 Set 

«(*) = / i t^tz 4 = / Y (t x 2 y 2 )(i + ^ + y 4 )- 1 dxdy. 

J x 2+y2<t l + x +y Jri 
First by Algorithm [TJ we compute a holonomic system for the integrand 

u(x, y, t) = Y(t ~x 2 - y 2 )(l + x 4 + y 4 )' 1 . 
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Then the integration algorithm outputs a holonomic system 
((t 5 + 3f 3 + 2t)d 2 + (2t 4 + 3t 2 )d t )v(t) = 0. 
Solving this differential equation by quadratures, one gets 



V (t) = & f - A = + c 2 

Jo Vs 4 + 3s 2 + 2 



Vs 4 + 3s 2 + 2 

for t > with some constants C\ , C^. On the other hand, by the change of 
variables to polar coordinates, one has 

v(0) = 0, lim v'(t) = ir. 

In conclusion, we have proved the identity 

ft ds 



v(t) = V2nY(t) f -= 
Jo v s 



dxdy. 



+ 3s 2 + 2 

Example 8 Set 

v(t)= [ dxdy= [ Y(t-x 6 -x 4 y 2 - y A ) 

Jx 6 +x 4 y 2 +y 4 <t JR 2 

By integrating a holonomic system for the integrand, we obtain a differential 
equation 

((147456t 7 - 995328t 6 )<9j + (3096576t 6 - 15925248t 5 )<9 t 6 

+ (20604416t 5 - 74822400t 4 )^ + (51215360t 4 - 115430400i 3 )9 t 4 
+ (43401540t 3 - 4677096(tt 2 )<9 t 3 + (8707020t 2 - 2078400t)d t 2 
+ (110880* - I05)d t )v(t) = 0. 

This differential equation has a regular singularity at t = with exponents 0, 
T5> 13' T3>T5> !§'!§• ^ seems difficult to identify v(t) in the 7-dimensional 
solution space of the above equation. 

Example 9 The integral 



v(t) = I e ~t(x 2 +y 2 ) dxdy = f Y(x 3 - y 2 )e- t{x2+y ^ dxdy 

Jx 3 -y 2 >0 JR 2 

satisfies a differential equation 

(216t 4 9 t 4 +(32i 4 +1836t 3 )a t 3 +(224t 3 +3594i 2 )9 2 +(326i 2 +1371t)a t +70t+15)w(t) = 
for t > 0. 



5 Difference-differential equations for definite in- 
tegrals 

Let us consider the case where the integrand has some auxiliary parameters 
a = (ai, . . . ,a p ) other than x = (xi,--- ,x n ). The general algorithm given 
in the preceding section does not work directly since the integrand does not 
necessarily satisfies a holonomic system including the parameters a as variables. 
To avoid this drawback, we consider difference-differential equations including 
the parameters. 
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Definition 5 A subset ft of C p is said to be shift-invariant if a e ft implies 
that a + (1, 0, . . . , 0), . . . , a + (0, . . . , 0, 1) also belong to ft. 

Let u(x, a) be a distribution in x — (x\ , . . . , x n ) on an open set U of W 1 with 
parameters a — [a\, . . . ,a p ) which belong to a shift-invariant subset ft of C p . 
Then the shift operator E ai with respect to aj acts on u(x,a) by 

E ai u(x, ai, . . . ,ai-i,a,i, a i+ \, . ..,a p ) = u(x, a\, . . . , aj_i, a^+1, o i+ i, . . . ,a p ) (i 

Let D„ be the ring of differential operators on the variables x. We denote by 
D n (a,E a ) the ring of difference-differential operators which is generated by ai 
and E ai (i — 1, . . . ,p) over D„ with the commutation relations 

-Ea^i — a,iE ai = { Q - |j ^ j| ' E ai E aj = E aj E ai , aiaj = ajai, 

where we assume that aj and E ai commute with the elements of D n . We 
introduce new variables t — (ti,...,t p ) and the associated derivations d t — 
(d tl , ■ ■ ■ , dt p ) and consider the ring D n+p of differential operators on the vari- 
ables (x,t) = (x 1 ,...,x n ,t 1 ,...,t p ). Let n : D n+p -> D n (a, E a , E^ 1 ) be the 
homomorphism of D„-algebra defined by 

This homomorphism is well-defined since 

- Ud u ) = fi(d ti )fJ,(ti) - fj,(U)fj,(d ti ) = -aiE~^E ai - E ai {-ai)E~^ = 1. 

It is easy to see that /i is injective and can be extended to an isomorphism of 
D n+p [(t 1 ...t p )- 1 } to D n (a, Ea.E^ 1 ). 

The inverse shift E~ x 'acts' on u(x, a) by 

E~\{x,a\, . . .,ai-i,ai,a i+ i, ...,a p ) = u(x,ai, . . . ,Oj_i,Oi - l,a i+ i, ...,a p ) 

if (fli, Oi — 1, ... , a p ) e ft. In general, an element P of D n (a, E a , E^ 1 ) acts 
on u(x,a) if a — fc(l, . . . , 1) £ ft for a sufficiently large integer k. Then an 
element Q of D n+p acts on a) by Qu(x,a) := n{Q)u{x,a). 

In order to justify working in the extended ring D n (a, i? a , -E^ 1 ), or in D n+p , 
let us make the following 

Assumption: Let u(x, a) be a distribution in x with parameters a defined 
on V x R d x ft with an open set V of M" _<i and a shift-invariant set ft of C p . 
The integral 

v(x',a):= / u(x, a) dx n - d+ i ■ ■ ■ dx n 

with x' = (xi, . . . , x n -d) is well-defined for (x',a) £ V x Q. Let P be an 
arbitrary element of D n -d(a, E a , E^ 1 ) such that Pv(x',a) = if x' € V and 
a — k(l, . . . , 1) e ft with a sufficiently large k. Then for any a£f such that 
E%P e D n _ d (a, E a ), one has E%Pv{x', a) = for x'el^andaG ft. 

This assumption is satisfied if w(a;', a) is holomorphic in a which belongs to a 
shift-invariant open set ft by virtue of the uniqueness of analytic continuation, 
or else if ft = W . 
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Definition 6 A left ideal I of D n (a 7 E a ) is called a holonomic D n (a, E a ) -ideal 
if 

^~ 1 {D n {a, E a , E^I) 
is a holonomic ideal of D n+p . 

Lemma 3 If a left ideal J of D n+p is holonomic, then /j,(J) n D n (a,E a ) is a 
holonomic D n (a, E a ) -ideal. 

Proof: It is sufficient to show the inclusion 

(D n (a, E a , E' 1 ) {fJ,{J) fl D n (a,E a ))) D J. 

If P belongs to J, then there exists a positive integer k such that 

Q ■■= E k ai ■ ■ ■ E k ^(P) = »(t k ■ ■ ■ t k p P) € D n (a, E a ) n M (J). 

Hence (i(P) = E~ k ---E- p k Q belongs to D n (a, E a , E' 1 ) (^(J) n D n (a, E a )). 
This completes the proof. □ 

Now let us describe an algorithm to compute a holonomic difference-differential 
system for the integral v(x', a). 

Algorithm 5 (difference-differential equations for an integral) 

Input: A set Go of generators of a holonomic ideal J of D n+p annihilating 
u(x, a). 

Output: A set G of generators of a holonomic ideal of D n _j(a, E a ) 
annihilating the integral v(x',a). 

1. Compute a set G\ of generators of the integration ideal 

N := D n _ d (t,d t ) n (d n - d+1 D n+p H h d n D n+p + J) 

of J with respect to x n -d+i, ■ ■ ■ ,x n . 

2. Let P be an element of G\. Then there exists a minimal 

v = (v\, . . . , Up) € W such that Q := E^/i(P) belongs to D„_d(a, E a ). Let us 
denote this Q by nm(/i(P)). Set G := {nm(/x(P)) | P e G*i}. 

Theorem 2 The left ideal of D n _d{a, E a ) which is generated by G is holonomic 
and annihilates v(x',a) for (x 1 , a) £ V x O. 

Proof: If P belongs to G\ , then there exist Qi € D n+p and R € J such that 

n 

i—n — d-\-l 

There exists v e W such that Q := £>(P), E v a ii(d Xl Q l ) and all belong 

to D n (a, E a ). Then one has 



)v(x',a) = / ' 



Qv(x,a)= / Qu(x,a)dx n -d+i---dx n 



n 

= T2 d Xi E^fj,(Qi)u(x,a)dx n - d+ i---dx n = 



i— n— d+1 
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if a — k(l, . . . , 1) g fl for a sufficiently large k. Then it follows from Assumption 
that [i(P)v = holds for any a £ ft. The holonomicity of the output follows 
from Theorem [1] and Lemma [3] □ 

As a typical example, let u be a holonomic distribution defined on an open set 
U of R" and fi, ■ ■ ■ , f p G R[x] be polynomials such that the set {x € J7 | fi(x) > 
(i = 1, . . . is non-empty. We assume that the product u(fi) s _^ ■ ■ ■ (/ p )+ is 
well-defined as distribution on R™ if s belongs to a shift-invariant open set f2 of 
C p . This is the case if u is a locally integrable function. We set ti = E Si and 
<9 ti = —SiE~ for i = 1, . . . ,p in the sequel. 

Algorithm 6 (Computing differential equations for u(/i)+ • • • ) 

Input: A set Go of generators of a holonomic ideal I of D„ annihilating a 
distribution and polynomials fx, ■ ■ ■ , fp G R[x]. 

Output: A set G of generators of a holonomic ideal J of D n+p annihilating 

"(A)? •••(/„)?. 

1. For P = P(:z,a ci) ...,a c J gG , set 



T ( P):=P|x,& l+ ^ ftj> ... >flUn+ ^ ftl 

This substitution is well-defined in the ring D n+p since the operators which 
are substituted for d Xl , . . . , <9 Xll commute with each other. 
2. Set 

G := {r(P) | P G G } U {i, - fox) \j = l,...,p}. 
Now let us prove the correctness of this algorithm: 
Theorem 3 Let J be the left ideal of D n+p generated by the set 

{T(P)\Pzl}U{t j -f i (x)\j = \,...,p}. 

Then J is a holonomic ideal of D n+p and annihilates u(fx)+ • • ■ (/ P )+ ■ 

Proof: First note that J is generated by the set 

{r(P)\PeG }U{t j -f j (x)\j = l,...,p} 

since r is a ring homomorphism. Assume that u(x) is G°°. Then in view of the 
equality 

we have 

r(p)( w (/ 1 ); 1 • • • = (p w )(/ 1 ); 1 ■ • • (/ P )J = o (6) 

if P G / and s — . . . , 1) G ^ for sufficiently large f G N. If u is locally 
integrable but not G°°, then ^ makes sense on {x G U \ fx(x) ■ ■ ■ f p (x) ^ 0} 
since (/i)+ ■ • • (/ P )+ is G°° in x there. We can verify that ((SJ) holds on U by 
using the same argument as the proof of Lemma 2.9 of [10] ■ It is easy to see that 
tj — fj(x) annihilates u(fi)+ ■ ■ ■ (/ P )+ • Hence J annihilates • • • (/ P )+ ■ 



22 



Let us show that D n+p /J is holonomic. Since D n /I is holonomic, its char- 
acteristic variety Chav(D n /I) is an n-dimensional algebraic set of C 2n . By the 
definition, we have 

Chax(D n+p /J) 

9fj_ , fa/; 



C {(x, t, £, r) G C 2 ^' | a(P) (x, a+E^"--^« + E = (VP € J), 

/,(x) ;./ 1 ?)} 



/,(x) :./ 1 ?)}• 



Since the set on the last line is in one-to-one correspondence with the set 
Char(_D„/7) x C p , the dimension of Char(_D„ +p / J) is n+p, which implies that 
Dn+p/J is a holonomic module. □ 

Example 10 As one of the simplest examples, let us consider the gamma func- 
tion 

/■oo poo 

v(s)=T(s + l)= e~ x x s dx= / dx. 

JO J-oo 

Introducing t — E s and d t — —sEj 1 , we compute a holonomic system 

(d x +d t + l)u= (t-x)u = Q 

for u(x, s) := e~ x x s + . Then integrating with respect to x yields a holonomic 
system 

(dt + l)v(a) = 
for v (s), which can be rewritten as 

(E s - (s + l))v(s) = 0. 

Example 11 Set 



v(a) = [ e~ x2 - y2 (x 3 - y 2 ) s dxdy = f e^"^ (x 3 - y 2 ) s + 

Jx 3 - y 2 >o Jm. 2 



dxdy, 



which is well-defined for s belonging to the set 

SI := {s G C | x ? -v, ~\-v,-\-v (v = 1, 2, 3, . . . )}. 



The holonomic ideal generated by 

-x 3 + y 2 + t, d x + 3d t x 2 + 2x, d y + (-2d t + 2)y 

annihilates the integrand u(x,s) :— e~ x2 ~ y2 (x 3 — y 2 ) s + . The integral ideal of 
this ideal is generated by 

108f 2 9 t 5 +(-216t 2 +648t)a t 4 +(108t 2 -972t+627)a t 3 +(356t-606)a t 2 +(-64t+108)a t +32t-48. 
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Hence v(s) satisfies a holonomic difference equation 

(32£ s 4 + (64s + 16)£ s 3 + (-108s 3 + 32s 2 + 32s)E 2 
+ (-216s 4 + 540s 3 -390s 2 + 66s)£ s -3s(s-l)(s-2)(6s-5)(6s-13))v(s) = 

for s G Q. 

Differential equations with parameters can be extracted from the output of 
the integration algorithm as follows. 

Algorithm 7 (intersection with a subring) 

Input: A set Go of generators of a left ideal J of D n+p , the ring of differential 
operators on the variables (xi, . . . , x n , tx, . . . , t p ). 

Output: A set G of generators of the left ideal J n D n [s] of D n [s] under the 
identification Sj = —dtjtj for j = 1, . . .,p. 

1. Introducing new variables Uj, Vj for j = 1, . . . ,p, let h(P) € D n+p [u] be the 
multi-homogenization of P £ D n+p ; i.e., h{P) is homogeneous with respect to 
the weight — 1 for t j and Uj , and 1 for 9 tj . , for each j . 

2. Let N be the left ideal of D n+p [u,v] generated by the set 

{h(P) | P e G } U {1 - UjVj \j = l,... , P }. 

3. Compute a set G\ of generators of the ideal N n D n+p by eliminating u, v 
via an appropriate Grobner base. 

4. Since each element P of Gi is multi-homogeneous without u, v, there exist a 
monomial S in t,dt and an operator Q(s) £ D n [s] such that 

5P = Q(-a tl < l7 ...,~9 tp t p ). 

Let G be the set of such Q for each P G Gi. 

The correctness of this algorithm was proved as Proposition 4.3 in [TT)] . 

Example 12 Consider the definite integral 

v(z,s)= / (x + t) s dxdy= I Y(l - x 2 - y 2 )(x + zf + dxdy 

JD(t) JR 2 

with D{z) — {(x,y) el 2 | x 2 + y 2 < 1, x + 2 > 0}, which is well-defined at 
least for Re s > and for any zei We first compute a differential holonomic 
system for the integrand in the variables [x, y, z, t) with t — E s and dt = —sE^ 1 . 
Then the integration algorithm outputs the ideal generated by 

d z + d u {z 2 -2tz + t 2 - l)d t + s-t. 

By computing the intersection with the ring D z [s] with the identification s = 
—dtt, we get a linear differential equation 

((1 - z 2 )d 2 + (2s + l)zd z - s(s + 2))v(z, s) = 0. 
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6 A shortcut for tensor product computation 



Let us consider again a definite integral of the form 

v(xu---,x„-d) ■= / u(x)Y(f 1 )---Y(f m )dx n - d+1 ---dx r , 



for a holonomic function or distribution u(x). In Section 4, we gave an algorithm 
to compute a holonomic ideal annihilating the integrand by a tensor product 
computation. Since the tensor product computation can often be a bottleneck 
in the actual computation, let us describe an alternative method using the 
technique introduced in the preceding section. 

Algorithm 8 (Computing differential equations for ^(/i)^ 1 • • ■ (/ p )+ p ) 

Input: A set Go of generators of a holonomic ideal I of D n annihilating a 
distribution u(x), polynomials f\, . . . , f p £ R[x], and complex numbers 
Ai, . . . , Xp. 

Output: A set G of generators of a holonomic ideal of D n annihilating 

uc/ifr ■••(/„)}.'. 

1. Let G\ be the output of Algorithm [5] with input Go and {fi, . . . , f p }. Let J 
be the ideal of D n+p generated by G\. 

2. Compute a set Gi of generators of Jn D n [s] with the identification 
Sj = —dtjtj for j = 1, ... ,p by using Algorithm [7] 

3. Set G := {P(A) | P(s) G G 2 }. 



Theorem 4 TTie output of the algorithm above generates a holonomic ideal and 
annihilates u(fi)^} ■ ■ ■ (/ p )^ p ifu(fi)^ ■ ■ ■ is well-defined and holomorphic 

with respect to s on a neighborhood of s = A. 



Proof: It follows from Theorem [3] and the uniqueness of analytic continuation 
that G annihilates u(Zi)^ 1 • • • (/ P )+ p - 

Set M :— D n /I. Let us show that J coincides with the annihilator of 
[1] ® / s in M ®c[x] D n+p f s , where [1] denotes the residue class of 1 G D n in M 
and / s = (/i) Sl • • • (f P ) Sp - Note that M ®c[x] D n+p f s has a natural structure 
of left Z?„ +p -module. 

First, let P be in J. Then there exist Q\, . . . , Qi, . . . , <3;+ p € D n+p and 
P±, . . . , Pi € I such that 

/ p 

p = E Q< r ( p + E Q'+i • - /,•(»))• 

It follows that 

i p 

i=i j=i 
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Conversely, suppose that P £ D n+p annihilates [1] ® f s . We can rewrite P 
in the form 

aeN™,i/GNP j=l j=l 

P 

3=1 

with p a ^(x) e C[x] and Qj € £>„ +p . Setting P„ := EaeN»Po 1 "( !i; )^'i we S et 

o = p([i] ® r) = [iv] ® (sr/ s ) e m ® c[x] £> n+p / s . 

It follows that each P„ belongs to / since {d"f 3 } constitute a free basis of 
Dn+pf s over C[x] (See e.g., Lemma 6.11 of [15]). Hence we have 

p = J2 w»o + E • (<i - /*(*)) g ^ 

There is a natural homomorphism 

p : M ® c[x] D n [s}f s -> M Oc^j D n+P / S ~ D n+p / J 

of left £> n [s]-module which sends [1] <g) f s to the modulo class [1] in D n+p /J. 
Moreover we have 

p{M® C [x] D n [s]f s ) = D n [s]/{JnD n [s]) C D n+p /J. 

Thus p induces a homomorphism 

v 

p' : M® c[x] D n f x ^D n [s]/[{JnD n [s})+J2{sj-Xj)D n [s]J 

3=1 

of left D„-module, which sends [1]<8>/ A to the residue class of 1, where f x denotes 
the residue class of f s in Ajs]/7(ELi(sj - *j) D n[s]f s )- Let N be the left 
ideal of D n generated by the output of Algorithm [5J Then by the definition we 
have an isomorphism 

p 

D n /N ~ D n [s}/((JnD n [s}) + J2(s3 -^)D n {s}) 

3=1 

as D n -module. In conclusion, p can be identified with a surjective homomor- 
phism 

p' : M ® c[x] Z?„/ A ->■ £> n /iV. 
This implies that D n /N is holonomic in view of Proposition [1] □ 

Remark 1 The homomorphisms p or p' are not injective in general; for exam- 
ple, set / = X1X2 and M := D n /I with I being the left ideal of D2 generated 
by xi and 82- Then 

9 2 ([i] ® f s ) = [i] ® d 2 f s = -[1] ® x!a f / s = ® = 
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holds in M ®c[x\ D2(t,d t )f s but not in M <8>c[x] D 2 [s]f s . The annihilator of 
[1] ® f° in M <S>c[x] -D2/ is generated by x% and £2^2, while N for A = is 
generated by x\ and c?2. 

Example 13 The Bessel function J v (x) satisfies a holonomic difference-differential 
system 

{x 2 d 2 x + xd x +x 2 - v 2 )J v (x) = {x{El + 1 ) - 2(y + l)E v )J v {x) = (7) 

for x £ M and v £ C\{— 1, — 2,...} and holomorphic in za Let us consider the 
integral 



u(z,^,s) := 




y s J v (x) dxdy 



with -D(z) := {(x, y) e ]R 2 | x, y > 0, a; 2 + y 2 < z). It is easy to verify that 
Y[x)J v {x) also satisfies ([7]). First we compute a holonomic system for 

u(x,y,z,v,s) := Y(z - x 2 - y 2 )y s + Y{x)J v {x) 

in the ring of differential operators on the variables x, y, z, ti, t<i with t\ = E v 
and t 2 = E s by using Algorithm [5] Tensor product computation fails for this 
example because of complexity. Then by the integration algorithm, we get 
a holonomic system for v(z, v, s) in the ring of differential operators on the 
variables z,i>, s. Finally computing the intersection with the subring D z [v, s] we 
get 

(8z 2 d 3 z + (-8s + 8)zd 2 + (2z - 2v 2 + 2s 2 )d z -s- l)v{z, is, s) = 0, 

which has a regular singularity at t = with exponents I, (±v + s)/2 + 1. This 
differential equation is valid, at least, for z £ K, Re s > —1, and v 7^ —1, —2, 

Finally we present timing data (in seconds) of the computation of difference 
equations in v for integrals involving J v {x) by using Risa/Asir running on a 
computer with 3.06 GHz Core 2 Duo Processor and 8 Gbyte memory. Step 
1 refers to the computation of a holonomic system for the integrand by using 
Algorithm [S] Step 2 consists of the integration algorithm (Algorithm [3]). Tensor 
product algorithm does not stop in a reasonable time period for most of the 
examples in the table. 



integral 


Step 1 


Step 2 


total 


/ Y(y)Y(l -x- y)Y(x)J v (x) dxdy 


0.08 


0.03 


0.11 


j Y{y)Y{\ -x 2 - y 2 )Y{x)J v {x) dxdy 


0.2 


2.2 


2.4 


J Y(y)Y(l -x 4 - y 4 )Y{x)J„{x) dxdy 


3.4 


346 


349 


j Y(y)Y(z)Y(l -x-y- z)Y(x)J u (x) dxdydz 


0.27 


0.13 


0.4 


J Y(y)Y{z)Y{\ - x 2 - y 2 - z 2 )Y{x)J lJ (x) dxdydz 


0.13 


2.4 


2.5 
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